home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 138_01 / rk4.c < prev    next >
Text File  |  1985-03-09  |  3KB  |  118 lines
  1. /*  
  2.     Runge - Kutta order 4 Algorithm
  3.  
  4.     Creation date:    31-Oct-83
  5.     Author:        Mike Roberts
  6.  
  7.     Copyright 1983, 1984 (c) California Institute of Technology.
  8.     All rights Reserved.  This program may be freely distributed
  9.     for all non-commercial purposes but may not be sold.
  10.  
  11.     This algorithm is described in detail on page 205 of
  12.     Burden, Richard L.:  Numerical Analysis.
  13.  
  14.     To approximate the solution of the initial value problem
  15.         y'=f(t,y),  a<=t<=b, y(a)=alpha,
  16.     at (N+1) equally spaced numbers in the interval [a,b]:
  17.     INPUT endpoints a,b; integerg N; initial condition alpha.
  18.     OUTPUT approximation w to y at the (N+1) values of t.
  19.  
  20. Step 1:
  21.     Set    h=(b-a)/N;
  22.         t=a;
  23.         w=alpha;
  24.     Output (t,w).
  25. Step 2:
  26.     For i=1,2,...,N do Steps 3-5:
  27.  
  28.     Step 3:
  29.         Set    K1=hf(t,w);
  30.             K2=hf(t+h/2,w+K1/2);
  31.             K3=hf(t+h/2,w+k2/2);
  32.             K4=hf(t+h,w+K3).
  33.     Step 4:
  34.         Set w=w+(K1+2K2+2K3+K4)/6;    (Compute w[i].)
  35.             t=a+ih.            (Compute t[i].)
  36.     Step 5:
  37.         Output (t,w).
  38.  
  39. Step 6:
  40.     Stop.
  41.  
  42.  
  43. */
  44.  
  45.  
  46. #define FALSE 0
  47. #define TRUE 1
  48.  
  49. rk4(function,a,b,n,alpha,t,w)
  50.  
  51. double (*function)();    /* function giving f(t,y) */
  52. double a;        /* beginning of interval */
  53. double b;        /* end of interval */
  54. int n;            /* number of steps in interval */
  55. double alpha;        /* initial condition for y */
  56. double t[];        /* array for returning T[i] values */
  57. double w[];        /* array for returning W[i] values */
  58.  
  59. {
  60.     register int i;    /* counter for integration steps */
  61.     double h;    /* stepsize */
  62.     double time;
  63.     double yapprox;    /* approximation for y value  */
  64.  
  65.     /*  STEP 1:  Initialization  */
  66.  
  67.     h = (b-a) / (double)n ;        /* Compute stepsize */
  68.     time = a;            /* Initialize time */
  69.     yapprox = alpha;        /* Start with the approximation
  70.                        equal to the initial value */
  71.     
  72.     for (i=0; i<n; i++)        /* Main integration loop */
  73.     {
  74.         if(i)    /* if not first time, call the integrator */
  75.             rk4_1(function,h,time,&yapprox);
  76.         
  77.             /* Pass the function pointer, the h, time, and
  78.                 yapprox values, and the pointers to
  79.                 the current positions in the T and W
  80.                 matrices */
  81.  
  82.         time = a + h*(double)i;  /* compute time */
  83.         t[i] = time;         /* also save it */
  84.         w[i] = yapprox;         /* store value for function */
  85.     }
  86. }
  87.  
  88. /*  This is the RK4 integrator portion.  It performs one step of the
  89.     integration, and is called on each step from the RK4 loop.
  90.     
  91.     function = pointer to function to integrate
  92.     h = stepsize
  93.     time = current time location
  94.     yapprox = current w (function approximation)
  95.  
  96. */
  97.  
  98. rk4_1(function,h,time,yapprox)
  99.  
  100. double (*function)();        /* Pointer to function to integrate */
  101. double h;            /* Step size */
  102. double time;            /* Current time step */
  103. double *yapprox;        /* Current approximation of function */
  104.  
  105. {
  106.     double k1, k2, k3, k4;    /* Temporary values in RK calculation */
  107.  
  108.     k1 = h * (*function)(time,*yapprox);    /* Evaluate first approx */
  109.     k2 = h * (*function)(time+h/2.0, *yapprox+k1/2.0);
  110.                         /* Evaluate second approx */
  111.     k3 = h * (*function)(time+h/2.0, *yapprox+k2/2.0);
  112.                         /* And the third */ 
  113.     k4 = h * (*function)(time+h, *yapprox+k3);
  114.                         /* And the last one */
  115.  
  116.     *yapprox += (k1 + 2.0*(k2 + k3) + k4)/6.0; /* new approx */
  117. }
  118.